WRIGHT-FISHER DIFFUSION WITH NEGATIVE MUTATION 

RATES 



SOUMIK PAL 

Abstract. We study a family of n-dimensional diffusions taking values in 
the unit simplex of vectors with non-negative coordinates that add up to one. 
These processes satisfy stochastic differential equations which are similar to the 
ones for the classical Wright-Fisher diffusions, except that the "mutation rates" 
are now nonpositive. This model, suggested by Aldous, appears in the study of 
a conjectured diffusion limit for a Markov chain on Cladograms. The striking 
feature of these models is that the boundary is not reflecting, and we kill the 
process once it hits the boundary. We derive the explicit exit distribution 
from the simplex, and probabilistic bounds on the exit time. We also prove 
that these processes can be viewed as a "stochastic time-reversal" of a Wright- 
Fisher process of increasing dimensions and conditioned at a random time. 
A key idea in our proofs is a skew-product construction using certain one- 
dimensional diffusions called Bessel-square processes of negative dimensions 
which have been recently introduced by Going- Jaeschke and Yor. 



1. Introduction 

An n-leaf Cladogram is an unrooted tree with n > 4 labeled leaves (vertices 
with degree one) and (n — 2) other unlabeled vertices (internal branchpoints) of 
degree three (see below for a figure). The number of edges in such a tree is exactly 
2n — 3. Sometimes they are also referred to as phylogenetic trees. Aldous in [3] 
proposes the following model of a reversible Markov chain on the space of all n-leaf 
Cladograms which consists of removing a random leaf (and its incident edge) and 
reattaching it to one of the remaining random edges. 

For a precise description we first define two operations on Cladograms. More 
details, with figures, can be found in [3]. 

(i) To remove a leaf i. The leaf i is attached by an edge e\ to a branchpoint b 
where two other edges e2 and es are incident. Delete edge e\ and branch- 
point 6, and then merge the two remaining edges e2 and e% into a single 
edge e. The resulting tree has 2n — 5 edges. 

(ii) To add a leaf to an edge /. Create a branchpoint b' which splits the edge 
/ into two edges f^fz and attach the leaf i to branchpoint b' via a new 
edge f\. This restores the number of leaves and edges to the tree. 

Let T n denote the finite collection of all n-leaf Cladograms. Write t' ~ t if 
t' 7^ t and t' can be obtained from t by following the two operations above for 
some choice of i and /. Thus a T n valued chain can be described by saying: 
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Figure 1. A 7-leaf Cladogram 



remove leaf i uniformly at random, then pick edge / at random and re-attach i 
to /. If we assume every edge to be of unit length, then it also involves resizing 
the edge length after every operation. In particular the transition matrix of this 
Markov chain is 



This leads to a symmetric, aperiodic, and irreducible finite state space Markov 
chain. Schweinsberg [16] proved that the relaxation time for this chain is 0(n 2 ), 
improving a previous result in [3]. 

In his webpage [2] Aldous asks the following question: what is an appropriate 
diffusion limit of this Markov chain? The invariant distribution for the Markov 
chain on n-leaf Cladograms is clearly the Uniform distribution. It is known, see Al- 
dous [1 , that the sequence of Uniform distributions on n-leaf Cladograms converge 
weakly to the law of the (Brownian) Continuum Random Tree (CRT). Hence, it is 
natural to look for an appropriate Markov process on the support of the CRT which 
can be thought of as a limit of the sequence of Markov chains described above. At 
this point it is important to understand that the support of the CRT consists of 
compact real trees with a measure describing the distribution of leaves. These trees 
are called continuum trees. For a formal definition of these concepts we refer the 
reader to the seminal work by Aldous in pQ. However, for an intuitive visualization, 
one should think of a typical continuum tree as a compact metric space on which 
branchpoints are dense and all edges are infinitesimally small. This implies that 
the Markov process that mimics the operation of removing and inserting a new leaf 
on a continuum tree should not jump, in other words, we can call it a diffusion. 

A detailed description of this diffusion on continuum trees is forthcoming in Pal 
[13] . In this article we consider several important features of this limiting diffusion 
that are of interest by themselves and provide bedrock for the followup construction. 

Consider the branchpoint b in the 7-leaf Cladogram t in Figure 1. It divides the 
collection of leaves naturally into three sets. Let X(t) = (Xi, X 2 , X 3 )(t) denote 
the vector of proportion of leaves in each set. The corresponding number of edges 
in these sets are (2nX\ — 1, 2nX2 — 1, 2nX^ — 1). For example, at time zero in our 
given tree, going clockwise from the right we have X(0) = (3/7, 2/7, 2/7). 

Let S n denote the unit simplex 




if t' 



if t' = t. 





Xi > for all i and 




NEGATIVE WRIGHT-FISHER 



3 



Some simple algebra will reveal that for any point x = (#1, #2, #3) in S3, given 
X(t) = x, the difference Xi(t') —X\(t) can only take values in { — l/n,0, 1/n} with 
corresponding probabilities 



2n(l — xi) — 2 2nxi — 1 

^ =Xl 2n -5 ' ^P-"^' ^ = (1 - Xl) ^3^- 



Thus 



e ( x l( to - x l( t) 1 x ( t) = x) = l 2xi ~t7 l] "-hV 1 - 



(2) E ((Xi(t>) - X x (t)) 2 I X(t) = x) = ^ 4nXl(1 ~ n ^ ) ~ Xl ~ 1 

« ^2xi(l - sci). 

If we take scaled limits, as n goes to infinity, of the first two conditional moments 
(the mixed moments can be similarly verified), it is intuitive (and follows by stan- 
dard tools) that as n goes to infinity, this Markov chain (run at n 2 /2 speed) will 
converge to a diffusion with a generator 

i,j=l i=l 

The generator written as above is similar to the generator for the well-known 
diffusion limit of the Wright-Fisher (WF) Markov chain models in population genet- 
ics. The WF model is one of the most popular models in population genetics. This 
is a multidimensional Markov chain which keeps track of the vector of proportions 
of certain genetic traits in a population of non-overlapping generations. A good 
source for an introduction to these models is Chapter 1 in the book by Durrett [6 . 
For computational purposes one often takes recourse to a diffusion approximation, 
which, in its standard form, leads to a family of diffusions parametrized by n "mu- 
tation rates" . The state space of the diffusion is given by S n and is parametrized 
by a vector (Si,..., S n ) of nonnegative entries. A weak solution of the WF diffusion 
with parameters S = (Si, . . . , S n ) solves the following stochastic differential equation 
for i = 1, 2, . . . , n: 

(4) dJi(t) = - (Si - 5 Ji(t)) dt + ^*i, j (J)dp j (t), $o = Yl 

j=i i=i 

Here j3 = (j3\, . . . , j3 n ) is a standard multidimensional Brownian motion and the 
diffusion matrix a is given by 

(5) °i,j( x ) = V%i {Hi = 3} - V x i x i) ' 1 - i'O - n - 

We define the Wright-Fisher diffusion with negative mutation rates to be a fam- 
ily of n-dimensional diffusions, parametrized by n nonnegative parameters S = 
(Si, . . . ,S n ), which is a weak solution of the following differential equation: 

^ n n 

(6) dfjLi(t) = --(Si- Som(t))dt + *y]aij(fj,)dl3j(i), S = Y]Sj- 

3=1 i=i 

The initial condition /x(0) is in the interior of S n and the process has a drift that 
pushes it outside the simplex. We will show later that the process is sure to hit 
the boundary of the simplex at which point we stop it. In the next section we will 
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explicitly construct a weak solution of §6y . The uniqueness in law of such a solution 
until it hits the boundary follows since the drift and the diffusion coefficients are 
smooth (hence, Lipschitz) inside the open unit simplex. The law of this process 
will then be denoted uniquely by NWF(<5i, . . . , S n ). 

Equivalently this process can be identified by its Markov generator. Expanding 
da' and using the fact that J27=i x i ~ 1> we S et 

(7) A n = \j2^ (H* = J) ~ *i) -itl&- ^7- 

i,j=l J i=l 

which identifies ^ as the generator for NWF(l/2, 1/2, 1/2). 

In this text we focus on properties of NWF models as a family of diffusions on 
the unit simplex and explore some of their properties that are important in the 
context of the Markov chain model on Cladograms. 

Part (1). We show that, just like Wright-Fisher diffusions (see [E]), the NWF 
processes can be recovered from a far simpler class of models, the Bessel- Square 
(BESQ) processes with negative dimensions. A comprehensive treatment of BESQ 
processes can be found in the book by Revuz & Yor [15 . This family of one 
dimensional diffusions is indexed by a single real parameter 6 (called the dimension) 
and are solutions of the stochastic differential equations 

(8) Z(t) = x + 2 / y/\Z(s)\dP(s) + 6t, x > 0, t>0, 

Jo 

where (3 is a one dimensional standard Brownian motion. We denote the law of 
this process by Q e x . It can be shown that the above SDE admits a unique strong 
solution until it hits the origin. The classical model only admits paramater 6 to 
be non- negative. However, an extension, introduced by Going- Jaeschke & Yor [7], 
allows the parameter 6 to be negative. It is important to note that Q e x is the 
diffusion limit of a Galton- Watson branching process with a \6\ rate of immigration 
(for 6 > 0) or emigration (for < 0). 

In Section |3] we show that the NWF(5i, . . . , S n ) law, starting at (xi, . . . , x n ), can 
be recovered via a stochastic time-change from a collection of n independent pro- 
cesses with laws Q~ 25 % i = 1, . . . ,n, and dividing each coordinate by the total sum. 
For the corresponding discrete models this is usually referred to as Poissonization. 

In this article we utilize this relationship to infer several properties about the 
NWF processes. For example, we prove that these diffusions, almost surely, hit 
the boundary of the simplex. We derive the explicit exit density supported on the 
union of the boundary walls in Theorem [9] 

Part (2). We also prove an interesting duality relationship between WF and 
NWF models. To describe the duality relationship we let the NWF continue in 
the lower dimensional simplex when any of the coordinates hit zero. Thus, every 
time a coordinate hits zero, the dimension of the process gets reduced by one, 
and ultimately the process is absorbed at the scalar one. Such a process can be 
obtained by running a WF model with appropriate parameters that initially starts 
with dimension one and value 1. At independent random times, the dimension 
of the process increases by one, and the newly added coordinate is initialized at 
zero. Finally we condition on the values of the process at a chosen random time. 
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The resulting process, backwards in time and suitably time- changed, is the original 
NWF model. 

Part (3). The time that the NWF process takes to exit the simplex is a crucial 
quantity due to a reason which we describe below. We keep our exposition mostly 
verbal without going into too much detail since the details require considerable 
formalism from the theory of continuum trees and will be discussed elsewhere. In 
[13] we show how Part (1) points towards a Poissonization of the entire Aldous 
Markov chain which is simpler for considering scaled limits. The Poissonized version 
of the Markov chain on n-leaf Cladograms stipulates: every existing leaf has an 
Exponential clock of rate 2 attached to it which determines the instances of their 
deaths, and every existing edge has an independent Exponential clock of rate 1 
attached to it at which point the edge is split and a new pair of vertices (one of 
which is a leaf) is introduced. It is an easy verification that the rates are consistent 
with the BESQ limit that we claimed in Part (1) above. Hence, one would expect 
that the limit of the Poissonized chains on continuum trees, normalized to give 
a leaf-mass measure one, and suitably time-changed would give the conjectured 
Aldous diffusion. This is the strategy followed in [13] . 

Now, the Poissonized chain has some beautiful and interesting structure. Please 
see pQ for the details about continuum trees that we use below. A continuum tree 
T comes with its associated (infinite) length measure (analogous to the Lebesgue 
measure) and a leaf-mass probability measure which describes how the leaves are 
distributed on it. We will denote the length measure by Leb(T) and the leaf- 
mass probability measure by /i(T). Suppose we sample n iid elements from /i(T) 
and draw the tree generated by them, it produces an n-leaf Cladogram with edge- 
lengths (or, a proper n-tree, according to [1 ). Thus, by using the fact that the 
continuum tree is compact, one can approximate a continuum tree by a sequence 
of n-leaf Cladograms. 

Now consider an n-leaf Cladogram for a very large n, and further consider m in- 
ternal branchpoints. For example, in Figure 1, we have three branchpoints {a, 6, c] 
in a 7-leaf Cladogram. These branchpoints generate a skeleton subtree of the orig- 
inal tree and partition the leaves as internal or external to the skeleton. The com- 
ponents of the vector of external leaf masses grow as independent continuous time, 
binary branching, Galton- Watson branching processes with a rate of branching / 
dying 2 and a rate of emigration 1. Note that this is consistent with the diffusion 
limit as BESQ with = — 1. As the Markov chain (Poissonized or not) proceeds, 
there comes a time when one of these external leaf masses gets exhausted. When 
this happens, one of the internal branchpoints becomes a leaf. The distribution of 
every coordinate of external leaf-masses at this exit time is derived in Part (2). 
Until this time, supported on the skeleton, new subtrees can grow and decay. We 
show in [13] , the dynamics of the sizes of these subtrees on the internal part can be 
modeled as the age process of a chronological splitting tree. Chronological splitting 
trees are special kind of biological trees where an individual lives up to a certain 
(possibly non-Exponential) lifetime and produces children at rate one during that 
lifetime. Her children behave in an identical manner with an independent and iden- 
tically distributed lifetime of their own. The age process refers to the point process 
of current ages of the existing members in the family. More details about splitting 
trees can be found in the article by Lambert [10] . 
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When one of the internal vertices gets exposed, the above dynamics breaks down 
and we need to find a slightly different set of internal vertices to proceed. Hence, 
it is important to derive estimates of the times at which this change happens. We 
provide quantitative bounds on the value of this stopping time under the special 
situation of symmetric choice of parameters which is the case at hand. 

The article is divided as follows. Our main tool in this analysis is to establish 
a relationship between NWF processes and Bessel-square processes of negative di- 
mensions, much in the spirit of Pal [12] . This has been done in Section [3] where we 
also establish Theorem [7| The relevant results about BESQ processes have been 
listed in Section [2] Most of these results are known and appropriate citations have 
been provided. Proofs of the rest have been done in the Appendix. Exact compu- 
tations of exit density from the simplex has been done in Section [4] Estimates of 
the exit time have been established in Section [U 



2. Some results about BESQ processes 

The Bessel-Square processes of negative dimensions — 0, where > 0, are one 
dimensional diffusions which are the unique strong solution of the SDE 

(9) X(t) = x-0t + 2 [ y/X(s)dP(s), t<T , 

Jo 

where T is the first hitting time of zero for the process X and x is a positive 
constant. The process is absorbed at zero. We will denote the law of this process 
Q~ e just as BESQ of a positive dimension will be denoted by Q e x . 

The following collection of results is important for us. All the proofs can be 
found in the article by Going- Jaeschke and Yor [7]. 

Lemma 1 (Time-reversal). For any 6 > — 2 and any x > 0, Q~ 6 (Tq < 00) = 1, 
while, for > 2, one has Q x (Tq < 00) = 0. 

Moreover the following equality holds in distribution: 

(10) (X(T -u\u< T ) = (Y(u), u<L x ), 

where Y has law Qq +9 and L x is the last hitting time of x for the process Y . 
In particular, 

(i) both L x and To are distributed as x/2G, where G is a Gamma random 
variable with parameter (0/2 + 1). 

(ii) The transition probabilities pf(x, y) for x, y > satisfy the identity 

The following results have been proved in the Appendix. 
Lemma 2. The scale function for Q~° , > ; is given by the function 

s(x) = x e/2+1 , x>0. 

Moreover, 

(i) The origin is an exit boundary for the diffusion and not an entry. 

(ii) The change of measure 

x -e/2-i Q -e (x(^/ 2+1 l(.)) 
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on the a-algebra generated by the process up to time t is consistent for 
various t and is the law of Q^ +e . Thus, we say Q^ +6> is Q~ e conditioned 
never to hit zero. 

The previous fact is the generalization of the well-known observation that Brown- 
ian motion, conditioned never to hit the origin, has the law of the three-dimensional 
Bessel process. 

Lemma 3. Let {Z(t), t > 0} denote a BESQ process of dimension for some 
> 2. Then 

1 f f du 1 , 

hm - — ; — — / — - = , tor all t > 0. 

e^o log(l/e) J e Z{u) 0-2' J 

3. Changing and reversing time 

Our objective in this section is to establish a time-reversal relationship between 
NWF and WF models. 

Theorem 4. Let . . . , z n and #i, . . . , 6 n be nonnegative constants. Let Z = 
(Zi, . . . , Z n ) be a vector ofn independent BESQ processes of dimensions — #i, . . . , — r 
respectively, starting from (zi, . . . , z n ). Let ( be the sum Ym=i 
Define 

T % = inf {t > : Zi(t) = 0} , r = A? =1 T<. 
Then, there is an n-dimensional diffusion \i, satisfying the SDE in (|6| for 
NWF{6\/2, . . . , n /2), for which the following equality holds: 

r tAr ds 

(11) Zi(t At) = C(t A r)/i, (4C t ) , 1 < i < n, C t = . 

Jo C{s) 

Thus, in particular, equation (|6| admits a weak solution for all non-negative pa- 
rameters . . . , S n ). 

Proof. The proof is almost identical to the case of WF model as shown in [I2j 
Proposition 11] with obvious modifications. For example, unlike the WF case, the 
time-change clock is no longer independent of the NWF process. We outline the 
basic steps below. 

We know from ^ that 

dZi(t A r) = -0id(t A r) + 2 V / Z^ i (t A r), i = 1, 2, . . . , n. 

Define O = EIU Let ^CO = Z */CW for t < T - Then b y It6 ' s rule we S et: 



(12) dVi(t A r) = -C- 1 [0i - OoVi] d(t A r) + y/V^l - V z )dM z (t), 

where 

(is) dM.it) = (H* = ft ~ yfiw) W A r )' 

and {Mi) (t) =4C t . 

Let {p n , u > 0} be the inverse of the increasing function 4(7^. Applying this 
time-change to the SDE for Vi in (12), we get 

(14) dfjLi(t) = ~ [0i - OoVi] dt + vV;(l-/^WW, 

where W$ is the Dambis-Dubins-Schwarz (DDS; see [15, p. 181]) Brownian motion 
associated with M { . This turns out to be the SDE for NWF{6 1 /2, . . . , n /2). □ 
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Let 0i, #2, . . • , n be nonnegative and z±, Z2, • • • , z n be positive constants. For 
i = 1, 2, . . . , n define independent random variables (Gi, . . . , G n ) where G^ is dis- 
tributed as Gamma(#i/2 + 1). Let 

(15) Rz = W z ' * = 1 > 2 >---,n. 

Also, independent of (Gi, . . . , G n ), let Yi, Y2, • • • , Y n be n independent BESQ pro- 
cesses of positive dimensions (4 + 0i), (4 + 02), . . . , (4 + n ) respectively, all of which 
are starting from zero. 

For any permutation tt of n labels, condition on the event 

(16) R ni > R n2 > ... > Rn n , and let R* = R n2 . 

We now construct the following n dimensional process (Xi, . . . , X n ): 

(17) Xi(t) = Y % ((t - R* + , t > 0. 

Notice that at time t = 0, every is at zero except the 7Tith. 

Let S(t) denote the total sum process Y^i=\Xi(t). Note that S(t) > for all 
t > with probability one. Define the process 

(18) G t := t>0. 

The process Gt is finite almost surely for every t (unfortunately, we cannot define 
R* — R ni precisely because Ct will be infinity; see Lemma [3]). Let A denote the 
inverse function of the continuous increasing function 4G. That is 

(19) A t = inf {u > : 4G n > t} , t > 0. 

Lemma 5. There is an n- dimensional diffusion v such that the following time- 
change relationship holds 

(20) Ui{t) = ^{A t ), or Xi(t) = S{t)vi(4Ct) t>0. 
The distribution of v is supported on the unit simplex 

S n = {xi > : xi + x 2 + . . . + x n = 1} . 

Conditional on the values of Gi, . . . , G n and the process S, the law of v can be 
described as below. 

Let 7r be any permutation ofn labels. On the event R ni > R* = R n2 > . . . > R^ n . 
Let V2 < ... <V n be defined by 

Ay. = R* — R ni , or, equivalently 4G j r*_ j r 7f . = Vi. 

Note that V 2 = 0. 

For i > 2 and Vi < t < Vi+i, the process v is zero on all coordinates except 
(7Ti, . . . , 7Ti). The process v{k\^ . . . , 7^), given the history of the process till time Vi 
( and the Gi ; s and S) is distributed as the classical Wright-Fisher diffusion starting 
from 

— (X 7Tl , . . . , X n . ) {Ay. ) = — (X 7Tl , . . . , X n . ) (R* — R 7Ti ) , 
and with parameters (j 7Tl , . . . , 7^ . ) where 

7i = i /2 + 2, j = l,2,...,n. 
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Proof. The Gamma random variables Gi, . . . , G n are independent of the BESQ pro- 
cess Yi, . . . , Y n . Thus, conditional on G?i, . . . , G n , the vector of processes (Xi, . . . , X n ) 
has the following description. For 

R* — Riii ^ t ^ R* — R-Ki+i i i ^ 2, 

all coordinates other than the 7Tith, 7T2th, . . ., 7r^th are zero. And, (X 7ri , . . . , X^), 
conditioned on the past, are independent BESQ processes of dimensions (4 + 
6 ni , . . . , 4 + 6 n . ) and starting from (X ni , . . . , X ni ) (R* — R ni ) . 

Thus, on this interval of time, the existence of the process identifying its law 
as the WF law, and the claimed independence from the process 5, all follow from 
[12j Prop. 11]. The proof of the Lemma now follows by combining the argument 
over the distinct intervals. □ 

Lemma 6. Consider the set-up in fll5| ), ( [l7| ) ; and ( fl9| ). Let Z2, . . . , Z n be n 

stochastic processes defined such that {Zi(t), < t < R*} is the time-reversal of 
the process {Xi(t), < t < R*}, conditioned on Xi(R*) = Z{. That is, conditioned 
on Xi(R*) = Z{ for every i, 

Zi(t)=Xi(R* -t)=Y i (R i -t) + , for 0<t<R*. 

Then (Z\, . . . , Z n ) are independent BESQ processes of dimensions — #1, . . . , — n , 
starting from zi,...,z n , and absorbed at the origin. 

Proof. It suffices to prove the following: 

Claim. Let {Y(t), t > 0} denote a BESQ process of dimension (4 + 0) starting 
from 0. Fix a z > 0. Let T be distributed as z/2G, where G is a Gamma random 
variable with parameter (0/2 + 1). Then, conditioned on T = I and Y(l) = z, the 
time-reversed process {Y((l — s) + ),0 < s < 00} is distributed as Q~ 6 ', absorbed at 
the origin, conditioned on To = L Here To is the hitting time of the origin for Q~ e . 

Once we prove this claim, the lemma follows since the law of To is exactly z/2G. 
See Lemma Q] 

Proof of Claim. For the case of 6 = 0, this is proved in [I4j p. 447]. The general 
proof is exactly similar and we outline just the steps and give references within [14] 
for the details. 

For any e R, t > 0, x, y > 0, let Q° x ^ y denote the law of the BESQ bridge of 
dimension 0, length t, from points x to y. That is to say, if Y follows Q e x , then Q e x ^ y 
is the law of the process {Y(s), < s < t} conditioned on the event {Y(t) = y}. 

Now, BESQ bridges satisfies time-reversal [TH p. 446]. Thus, if we define P to 
be the P-distribution of a process {X(t — s), < s < t}, then Q e x ^ y = Qy'l+x- 

We consider the case when the dimension is (4 + 9), 6 > 0, x = 0, y = z > 0. 
Then 

n 4+6>,£ _ ^4+6>,£ 

Now, from Lemma |5] (also see [I4j Section 3, p. 440]) we know that Q^ +e is 
Q~° conditioned never to hit zero (or equivalently, Q~ e can be interpreted as Q^ +e 
conditioned to hit zero). Since the origin is an exit distribution for Q~ e and not an 
entry (Lemma [2] see [14J p. 441] for the details of these definitions), the conditional 
law Q^lfo* ^ s n °thing but Q~ 6 , conditioned on Tq = t. This completes the proof. □ 
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The following is a more precise statement. 

Let (zi, . . . , z n ) be a point in the n-dimensional unit simplex S n . Fix n nonneg- 
ative parameters . . . , 5 n . Let Gi, . . . , G n denote n independent Gamma random 
variables with parameters 5i + 1, . . . , S n + 1 respectively. Define Ri = zi/2Gi. 

For any permutation tt of n labels, condition on the event R ni > R n2 > . . . > 
R 7Vn , and let R* = R n2 . 

Define the continuous process S by prescribing S(0) = Zi(R ni — R*) where Z\ 
is distributed as Q 7ri ; and for any t such that 

i?* — R ni <t<R*— ii 7ri _ | _ 1 , i > 2, Rir n+1 — 0. 

given the history, the process is distributed as a Bessel-square process of dimension 
Ei=i( 4 + 2 M starting from 5(i?* - i^J. 
Define the stochastic clocks 

du — diz 



C t = / ^r, C t = / — , 0<f <iT. 

7 S(u) J R *_ t S(u) 

and let A t denote the inverse function of 4Cf Let V2 < . . . < V n be defined by 
ACr^-r^. = V{. Note that V2 = 0. The 4 is a standardization constant that appears 
due to the factor of 2 in the diffusion coefficient in (J8|. 

Define an n-dimensional process given . . . , R n and the process S. For i > 2 
and < t < Vi+i, the process ^ is zero on all coordinates except possibly at indices 
(7Ti, . . . , 7Ti). At time zero, the process starts at the vector that is 1 in the 7Tith 
coordinate and zero elsewhere. 

Conditioned on the history till time V*, the process {z/(7Ti, . . . , 7r^)(£), V* < 
£ < ^+1} is distributed as the classical Wright-Fisher diffusion starting from 
^(7Ti, . . . , 7Ti)(Vi) and with parameters (j 7Tl , . . . , j ni ) where 

7j = $3 + 2 > j = 1,2, ...,n. 

Finally, consider the conditional law of the process, conditioned on the event 

S(R*)vi(AC R *) = z u for all i = l,2,...,n. 

Theorem 7. Define the time- reversed process 

where o denotes composition. Then this conditional stochastic time-reversed pro- 
cess, until the first time any of the coordinates hit zero, has a marginal distribution 
(when Gi 's and S are integrated out) NWF(5i, . . . , 5 n ) starting from (zi, . . . , z n ). 

Proof of Theorem^ We start with given values of R ni > R n2 > . . . > R 7Vn and the 
process S. Applying equation (20) in Lemma[5]to obtain the processes (Xi, . . . , X n ) 
defined by 

Xi(t) = S(t)vi(4C t ), 0<t<R*. 
Then, the vector (Xi,X2, • • • ,X n ) has the law prescribed by ([l7|). 



Now we apply Lemma [6] to obtain (Z\, . . . , Z n ) by conditioning (Xi, . . . , X n ) 
and reversing time. Finally the construction in Theorem [4] gives us the vector 
(/ii, . . . , fi n ) from (Zi, . . . , Z n ) as desired. □ 
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4. Exit density 

Let Zi, Z2, . . . , Z n be independent BESQ processes of dimensions — #i, . . . , — 9 ni 
where each Oi > 0. We assume that at time zero, the vector Z = (Zi,...,Z n ) 
starts from a point z = (21, ... , z n ) where every ^ > 0. Define Ti to be the first 
hitting time of zero for the process and let r = A^T^ denotes the first time any 
coordinate hits zero. We would like to determine the joint distribution of (r, Z(r)). 

Note that since each Ti is a continuous random variable, the minimum is attained 
at a unique i. Thus, for a fixed 1 < i < n, conditioned on the event r = T^, the 
distribution of Zi(r) is the unit mass at zero, and the distribution of every other 
Zj(t) is supported on (0, 00). Now, let hi denote the density of the stopping time 
Ti on (0, 00) and let refers to the transition density of Q~° . It follows that for 
any aj > 0, j ^ z, we get 

P(r = Ti,r < t, Zj(r) > for all j ^ i) 

=P < t, Tj > Ti, Zj(Ti) > a,, for all j + i) 

— I hi(s) TT P (Tj > s, Zj(s) > <2j) ds = / hi(s) TT P (Zj(s) > aj) ds, since aj > 0. 
Jo • / • Jo • / • 

ds. 



: / hi(s) Y[ / q-^(z j ,y j )dy. 

JO Ja 3 



Our first job is to find closed form expressions of the integral above. To do this 
we start by noting that Ti is distributed as Zi/2Gi (see Lemma [T]), where Gi is a 
Gamma random variable with parameter (4 + 0i)/2 — 1 = 0^/2 + 1. That is, the 
density of Gi is supported on (0, 00) and is given by 

y e ' /2 -y 



IW2 + 1) 
It follows that 

h is) ~ (**/ 2 )* i/2+1 -flj/2-2 - gi /2s < 8 < OO 

A) T(0i/2 + l) S 6 ' u ^ 5<0 °- 

On the other hand, it follows from time reversal (Lemma hi) that q s dj (zj, yj) = 
q s 3 (yj, Zj). For any positive a, the transition density q^ (y, z) is explicitly known 
(see, e.g., in [12]) to be s~ x f(z/ s,a,y / s), where /(-,fc, A) is the density of a non- 
central Chi-square distribution with /c-degrees of freedom and a noncentrality pa- 
rameter value A. In particular, it can be written as a Poisson mixture of central 
Chi-square (or, Gamma) densities. Thus, we have the following expansion 

(21) =^(^^ 

k=0 

where g r is the Gamma density with parameters (r/2, 1/2). That is 



Now, define 



o-r/2 r/2-l 

*■(*> = -r^r e " x/2 ' x -°- 

j¥* j¥* 
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Thus 



m/2 + l) 8 f 11 * ^ fc! r(^/2 + 2 + fe) 



-di/2-2-(n-l) -Zi/2s \ 



(gi/2)W 

r(^/2 + 1) ' 

-(i7 4 +?i)/2s 2 -e i /2-2(n-i) n V 2- fc (z,/ g )^/ 2 + fc + 1 



We now exchange the product and the sum in the above. We will need some 
more notations for a compact representation. For any two vectors a and 6, denote 
by 



a h = \[a\\ a\ = Y[a t l 



Also let ®i,yi,Zi stand for the vectors (9j, j ^ z), j ^ z), and (zj, j ^ z) 
respectively. 

Let k denote the vector (fcj, j ^ z), where every fcj takes any nonnegative integer 
values. Let k'l be the sum of the coordinates of k. Then 

°° (y j /2s) k 2- k (z j /s)W 2+k + 1 



r((9 j /2 + 2 + fe) 



k „k 



~k ( } 1 k ^ k! n^mv2+2 +%) - 

Thus, combining the expressions we get 

0i/2+l 

M.)n«^te.»)-r^Tij"-* /a - 1 - v ^- 

(22) J ^ 

V 7 CO 

s -^/2-2-(n-l) e -2 ;i /2 Se -(^+^)/2s ^ ^-N s -Si/2-2N-(n-l) Bn ^ 

N=0 

where 

©i/2+1 



k'l=7V 11 



^r^/2 + 2 + fc,-)' 



We can now integrate over s in ( 22 ) to obtain 

/•CO 00 /'CO 

/ h i (s)\[ q - 9 >(z 3 ,y 3 )ds=Y J B' N s- a »e- b / s d S , 
Jo jfti N=0 Jo 
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where 




(23) 


b 'n 


(24) 


«jv 


(25) 


b 



0</2+l n 

m/2+i) 2 4 Bjv ' 

6»i/2 + 0^/2 + 2n + 27V = 6> /2 + 2n + 27V, 

n 

^/2 + + £)/2 = (j7< + * )/2, z = Xl^- 

i=l 

Now a simple change of variable w = 1/s shows 

rOO rOO nOO 

/ s -a Ne -b/s ds= I w a Ne -bw w -2 dw = / ^iv^g-ta^ 

iW-i) f° ^ ^"-*e-»°L - r{aN - x) 



ft^- 1 y r(aiv-i) ^ a7V_1 

Since the zth coordinate of the exit point is zero, one can define yi = and 
Vo — Y^j=i Vj = Vi t° simplify notation. Thus we obtain 

+ o -0n/2-2n+l,i-iV d r ( a ^ _ 1) 



0i/2+l ©i/2+1 oo 

_ g» Z i o-e /2-2»+l //-• , w 9 \-e /2-2n-2iV+l 

- m/2 + 1) A,(^ + ^/2) 

x r( V2 + 2n + 2. - 1)4- g ^ ^ + 2 + fcj ) 

^0/2+1 00 

_ 9 -6> /2-2n+l / _i \-6> /2-2n-2AT+l 9 6> /2+2n+2AT-l 

"W2 + 1) 2 Ifo ° } 



x F( V2 + 2n + 2^ - 1)4- g n^fel 



) 



_©/2+l 00 

z / x-6> /2-2n-2iV+l 



m/2+i) ^ 



r(V2 + 2n + 2^-l) k E g ^ r(gj f 2 + 2 + fcj) 



k'l=iV 

We have the following result. 

Theorem 8. Let Zi,Z2,...,Z n be independent BESQ processes of dimensions 
—0i, . . . , — n , where each Oi > 0. Assume that Zi(0) = 2^(0) > 0, for every i. 

The distribution of (r, Z(r)) zs supported on the set (0, oo) x Uf =1 Hi, where Hi 
is the subspace orthogonal to the ith canonical basis vector e^. T7m£ zs 

#i = {(2/1,2/2, — , • 2/i = 0}. 

(i) Let Gi, i = 1, 2, . . . n be independent Gamma random variables with param- 
eters 0^/2 + 1, i = 1, 2, . . . , n. The law of t is the same as that o/min^ ^fjr 
and 

P( T = T i )=p(^> **>-, /or all j ± i) , 
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(26) 



where Ti is the first hitting time of Hi . 
(ii) The restriction of the law of the random vector Z{r), restricted to the hy- 
perplane Hi, admits a density with respect to all the variables y 3 - 's, j ^ i, 
which is given by 

Ql-6> /2-2n 71 n , 00 

= ^TI) II E r«V2 + 2n + 2N- 1)5-- 



k.AT(Q. 

E, V ; kj=Nj^i 



E n %! r(^./2 + 2 + ^)' 



i7ere 

n n 

5 = + m = °' °° = XX 

2 = 1 2=1 

Using Theorem [4j we get that the exit distribution of iVWF(5i, . . . , S n ) starting 
from a point (^i, . . . , z n ) G S n is the image under the map 

Xi H> =^ , 1 < i < n 

2^3 = 1 X 3 

of the exit density of independent BESQ processes of dimensions —Q\, . . . ,—Q n , 
where each 6i = 25 i. 

Theorem 9. The exit density of \i ~ NWF(5i, . . . , 5 n ) starting from {z\, . . . , z n ) G 
S n is supported on the set \Jf =1 Fi f where Fi is the face {x G S n : x % ; = 0} ; and 
admits the following description. 

(i) Let Gj, i = 1, 2, . . . n be independent Gamma random variables with param- 
eters 5i + 1, z = 1, 2, . . . , n. Then 

I G G 

(27) P (/i exifo through Fi) = P I — >—?-, for all j ^ i 

(ii) Let 5 represent the vector . . . , 5 n ) and Ze£ Jo = SlLi ^* ^ e 
tribution of the process \i, restricted to Fi, admits a density with respect to 
all the variables Xj 's, j ^ i, which is given by 

(28) (^ + l)E r r ( y+2n+t ) ) S ^„(z;k + ,5 + 2)^r„_ 1 ( a; ;k+l). 



AT=0 



i7ere the inner sum above is over all nonnegative integers (kj, j ^ i) such 
that ^Zj^i kj = N. The vector k represents a vector whose jth coordinate 
is kj for all j ^ i, and k 2 = 0. The vectors k + 5 + 2 and k + 1 represent 
vector additions of k, 5, and t/ie vector of all twos, and k and t/ie vector 
of all ones respectively. The factor Dir n -\ is a density with respect to 
the (n — 1)- dimensional vector (xj, j 7^ i) with corresponding parameters 
(kj + 1, j ^ i). It can also be interpreted as the conditional density of the 
n-dimensional Dir n (x; k + 1), conditioned on xi =0. 



Note that the density in ([28]) is a mixture of Dirichlet densities, strikingly sim- 
ilar to those appearing as transition probabilities of the Wright-Fisher diffusions 
themselves. See Griffiths [8 j, Barbour, Ethier, and Griffiths [4, and Pal [T2] . 
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Proof of Theorem [9j This is a straightforward integration. We have assumed that 
T,i z i = !• Tnus > S = 1 + J2j Vj\ define Vo = Ej 2/i> and 

x j=yj/yo, i<j<n. 



Hence (|26j) simplifies to 

2-2n 

II z i E r (^/ 2 + 2n + 2iV - 1)(1 + y ) 

(XjZj) 



(^ J_ U-0o/2-2n n „ , 

= (1 + 2/p) rr ^i/2+i V" Tva_ /o i o,„ i o/u 1 vi , ,,^-2V 

(29) 



x 7/ N V IT 



Now, to get to formula (28) we need to make a multivariate change of variables. 
Without loss of generality, let i = n. Then, for any y G Fi : we have i/ n = 0. Define 
the change of variables 

(2/1, . . . , 2/ n _2, 2/n-l) ^ (2/0, El, • • • , Xn-2) • 

In other words, yi = y^i, for alH = 1, 2, . . . , n— 2 and 2/ n _i = ?/o(l — ^i — • • .-^-2)- 
The determinant of the well-known Jacobian matrix is given by y^~ 2 • 
Thus, the density of . . . , x n ) restricted to Fi is given by 



(30) 

n 

w?iy n^ +1 Em/2 + 2n + 2iV-l) 

V / ; 7 j=l AT=0 



V +n " 2 (l + yf- d « l2 - 2n - 2N dy x ^ J] 



Jo * v— ^/ - l± k \y(0 ./2 + 2 + jb 7 -)" 



The following formula is easily verifiable for a > 0, /? > a + 1: 



Jo 



y a (l^y)- (3 dy = [ x^' 2 ^ - x) a dx = B(a + 1, f3 - a - 1) 



Jo Jo 

where 5 refers to the Beta function. 
In other words, (30) reduces to 



n 00 



7V=0 



(31) W + n-l,» + n + ,.m £ n w ^ t¥ 

T,&i*i=N&i 3X31 3J 

We now change 0^/2 to ^ and rewrite the above expression in terms of Dirichlet 
densities. We use the notations in the statement of Theorem |9l the vector k 
represents a vector whose jth coordinate is kj for all j 7^ z, and = 0. The 
vectors k + S + 2 and k + 1 represent vector additions of k, (5, and the vector of all 
twos, and k and the vector of all ones respectively. The factor Dir n _i is a density 
with respect to the (n — l)-dimensional vector (xj, j ^ i) with corresponding 
parameters (kj + 1, j ^ i). It can also be interpreted as the conditional density of 
the n-dimensional Dir n (x; k + 1), conditioned on ^ = 0. 
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Hence, for any (kj, j ^ i), integers, 

6i + l y k 3+ 5 3+ 1 ^3 

z i TT z 3 ZJ_ = 

^ + 1)11^ + 2 + ^) k \ 



{8t + _i) Dir n (z; k + J + 2)Dir n _ 1 (x; k + 1). 



r(5 + N + 2n)r(7V + n - 1) 
Thus (|3l|) reduces to 

r(6 + 2n + 2N -l)B(N + n-l,N + n + 6 ) 



oo 



(32) 



{Sl + 1} £ r(5 + + 2n)r(7V + n - 1) 

^ Dir n (z;k + (5 + 2)Dir n _i(x;k + l) 



x 

k'l=7V 

However 

r(^ + 2n + 27V - l)£(iV + n - 1, iV + + £ ) 
rpo + JV + 2n)r(7V + n - 1) 
r(J + 2n + 27V - 1) r(7V + n - l)T(N + n + S ) 
~~ T(5 + TV + 2n)r(7V + n - 1) r(27V + 2n + <5 - 1) 
= r(7V + n + go) 
~ r(7V + 2n + ^o)' 



This completes the proof of formula (28). 

The probability in (27) is a direct consequence of Theorem [8] conclusion (i). □ 

5. Exit time 

Let X = (Xi, . . .,X n ) be distributed as NWF(-0i/2, . . . , -6 n /2) starting from 
a point (xi, . . . , x n ) in the unit simplex. Let <jq denote the stopping time 

<Jo = inf {t > : Xi = for some i} . 

Our objective is to find estimates on the law of cfq. 

We will simplify the situation by assuming that all xi = 1/n and all Oi = 6. To 
this end we use the time-change relationship in Theorem[4j Let Z = (Z\, . . . , Z n ) be 
independent BESQ processes starting from (zi, . . . , z n ) as in the set-up of Theorem 
[4j where each z\ is now one. Then 

(33) a =4 / CO0 = £W 

Jo . =1 

By Theorem[8j the distribution of r is the same as considering n iid Gamma(#/2+ 
1) random variables G\, . . . , G n , and defining 

(34) T = 2^G? 

Our first step will be to prove a concentration estimate of max^ G^. 

Lemma 10. Let Gi, G2, . . . , G n be n iid Gamma random variables with parameter 
r/2, for some r > 2. Let x be the random variable max^Gi. Then, as n tends to 
infinity, 



E\[X = © ( V lo g n 
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Proof. First let r G N. Let {Zi(i), ...,Z n (i), i = 1, 2, . . . , r} be a collection of iid 
standard Normal random variables. Then 2Gj has the same law as Z 2 (l) + ... + 
Z?(r). Hence 

£max |Z|, (1) < E^x < Vr~Emax \Z\. (i) 

As n tends to infinity, the right hand side above converges to ^2r log(rn) while 
the left hand side converges to y/2 logn. This completes the argument for r G N. 
For a general positive r, bound on both sides by [r\ and |_ r J + 1- ^ 

We also need a version of logarithmic Sobolev inequality for Gamma random 
variables which can be found in several articles including [5]. 

Lemma 11. [5 , p. 2718] Let \i d denote the product probability measure of n iid 
Gamma(0) random variables. Then, for every f on W 1 which is in C 1 (i.e., once 
continuously differentiable) one has 

(35) Ent(/ 2 ) <aJ \jt Xi (d t f(x))^j d^(x). 

Here Ent(-) refers to the entropy defined by 

Ent(/ 2 ) = J / 2 log(/W " (/ /V) log (| fW 

And di refers to the partial derivative with respect to the ith coordinate. 



Lemma 12. Consider the set-up in Lemma [11} Let F be a function on the open 
positive quadrant (i.e., every Xi > 0) which is C 1 and satisfies 

n 

(36) (d;F) 2 < F. 

2 = 1 

Then the following concentration estimate holds for any r > 0: 

// (Vf - E 9 Vf > r) < exp (-r 2 ) , // (Vf - E 9 Vf < -rj < exp (-r 2 ) 
where E q VF = j VFdfi 6 . 

Proof. Condition ( [36] ) implies that 4 J^ =1 xi (di\[F^ < 1. Hence, from the classi- 
cal Herbst argument (for example, the monograph by Ledoux [H]), with a gradient 
defined by the right side of ( 35 ) we get 



(y/F-E e VF>rj < exp I 



-r 2 ) 



Here ^(y/F) is the expectation of y/F under \i e . Repeating the argument with 
— \/F instead of \[F we get the result. □ 

Theorem 13. The random variable x — max^G^ where G^s are iid Gamma(O), 
satisfies the following concentration estimate 

(37) P ( y/x > E( y/x) + r) < e~ r2 , for all r > 0. 
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Proof of Theorem [13} To prove (37) we start by noting that Lemma 12 is satisfied 
by the family of L^norms, {F kl k > 1}, defined by 

/ n \ iA 



F k (x) 



\i=i 



This is because each F k is smooth (when every Xi is positive) and 

2 



(38) ^.r ; uV, /•;,(•'•))- £\ 



1-1//C 



n 2/c-l 



(e; =1 4) 



2-2/fc ' 



Since for any nonnegative i/i, 2/2? • • • > 2/n an d an y > 1 one nas 

n / n \ P 



applying it for ^ = x J 5 and /3 = 2 — we get 

n / n 



2-1/fc 



i=l 



v?;=i 



Combining the above with ( |38| ) we get 

n / n 



i/k 



F k (x). 



Thus Ffe satisfies condition (36). 

Since F k converges pointwise to max^ xi as k tends to infinity, by applying DCT, 
Lemma 12 is true for the function max^ G^. This proves (37). □ 



Our next step will be to prove estimate on the quantity <Jo in ( 33 ) . The process 
£(s) is non-Markovian and not distributed as Q~ n6 '. However, on an possibly en- 
larged sample space one can create a Q~ n0 process ( such that the paths of ( and 
( are indistinguishable until ctq. This is possible by considering the SDE solved by 





((t) = n-n6t 



Jo 



t < a . 



To extend the process beyond 00, one concatenates an independent Brownian mo- 
tion W and defines 



(W(t), t < a 

\W(t) + W(t-a ), t > <tq. 



Then j3 is a Brownian motion in the enlarged filtration. Since Q n0 admits a strong 
solution, the process 



(39) 



C(t) =n-n0t + 2 



((s)dW(s), t < T , 



has law Q~ n0 and pathwise indistinguishable from £ until time do- Thus in the 
following discussion we will treat as if ( itself is distributed as Q~ n0 keeping in 
mind the above construction. 
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Theorem 14. Let n be distributed as an n- dimensional NWF(5, 5, ... ,5) starting 
from the point (1/n, 1/n, . . . , 1/n). Let ctq be the first time that any of the coordi- 
nates of \i hit zero. Let 



E max \J~G~i, Gi % ~ Gamma(S + 1). 

Ki<n 



Then, a n — 0(vTogn) ; cr has the law given by ( J33| ) where ( is distributed as Q\ 
and r is a random time. 

Moreover, for any r > 0, we get 

1 KV^K—-^ ) > 1 — 2e" 



2n5 



n(a n + r) - n(a n + r) 

Remark 1 . It is impossible to provide a simple description of the exact distribution 
of <Jo due to the distributional dependence of £ and r. The above theorem shows 
that r is about a constant and one can compare the distribution of o~o with that 
of 4 J Q du/((u) : where the upper limit of the integral is a constant. Limiting large 
deviation behavior of such integrals is possible to derive by methods as in [T7] . 



Proof of Theorem 14 , The proof is obvious from Lemma 13 and expression ( 34 ). □ 



Appendix: 

Proofs of properties of BESQ processes 

Proof of Lemma^ We use exercise 3.20 in [15j p. 311]. The scale function for 
Q e for > is well-known to be x _6> / 2+1 (see [15j p. 443]). Nearly identical 
calculations lead to the case when is replaced by —9 and we obtain the scale 
function s(x) = x 6>//2+1 . 

The speed measure is the measure with the density 

I* x 2 1 Q/2, 1 

m {X) = 7(xjAx = 2{9/2 + l) X 

We now Feller's criterion to check if the origin is an entry and/or exit point (see 
p. 108]). Note that 

(40, 1/2 > = whv f - to - 2,/a ) 

ra(0, £] = oo, for all positive £. 

Thus 

,1/2 ,1/2 

/ ra(f,l/2]s(df) < (X), and / ra(0,f]s(df) = oo. 

This proves that the origin is an exit and not an entry. 

Finally, to obtain part (ii) we apply Girsanov's theorem ([TSJ p. 327]). Let X 
satisfy the SDE dX(t) = -Odt + 2y/X(t)d/3(t), then, we take D(t) = X e ^ 2+1 (t) 
(without the normalization, for simplicity) and apply Girsanov. Under the changed 
measure, there is a standard Brownian motion /3* such that 

m = F'(t)+ f x- e ' 2 - i ( S )d{p,D) s 

Jo 

= P*(t)+ [ X-°/ 2 - 1 (s)(6 + 2)X e / 2+1 / 2 (s)ds = p*(t) + (Q + 2) I X-^ 2 (s)ds. 
Jo Jo 
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Thus, under the changed measure 

dX(t) = -Odt + 2X 1/2 (t)df3(t) = -Odt + 2(0 + 2)dt + d/3*(t) = (0 + 4)dt + d/3*(i). 

The interpretation as the conditional distribution is classical (see [H]). □ 

Proof of Lemma [3J For the assertion it is enough to take t = 1. Note that, under 
Qq, the coordinate process satisfies time-inversion, i.e., the process {t 2 Z(l/£), t > 
0} has law Qq. Thus, for < e < 1, if we define 

^ dt 



r 1 ju_ _ r 1 



t 2 z(i/ty 



then U e has the same law as Ci/ e — C\ = J^ e du/Z(u). Thus, by [T71 Thm. 1.1], 
we get lim e ^o U e / log(l/e) = (0 — 2) _1 almost surely. □ 
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